function Collect_t_pigou_Heterinc_m_pos(bs_max)

%bs_max = 16 

path_data = {'/Users/shoude/Dropbox/eegap/EEgap_data_code_heter_SM/Data/Matlab_estimation/'};
path_results = {'/Users/shoude/Dropbox/eegap/EEgap_data_code_heter_SM/Results_Simulation/'};


%Code Map
%For each income group we have evaluated the welfare for different values of the standard
%We have also done for each bs iteration.
%This code take different combination of bs iterations, and aggregate welfare by weighting each income group.

%   use $pathname\refrigerators\lcidemo_046_2008_2011_struct_v11_robustb, clear
% tab group_id if group_id>0 & group_id<7
% 
%    group_id |      Freq.     Percent        Cum.
% ------------+-----------------------------------
%           1 |    224,874       12.23       12.23
%           2 |    309,652       16.83       29.06
%           3 |    464,059       25.23       54.29
%           4 |    335,493       18.24       72.53
%           5 |    217,245       11.81       84.34
%           6 |    288,036       15.66      100.00
% ------------+-----------------------------------
%       Total |  1,839,359      100.00


tmp_file=strcat(path_data{1},'share_inc_heter.csv');
share_inc_heter=readtable(tmp_file);	    
%share_inc_heter(:,1)=[];  
share_inc_heter=table2array(share_inc_heter)'; 

bs_1_16=[1:16];
all_comb_bs_inc=nmultichoosek(1:16,6);
s = RandStream('mlfg6331_64');
bs_inc= datasample(s,all_comb_bs_inc,bs_max,1,'Replace',false);

%for i=1:150
%		 carbon_price(i)=2*(i-50);
		 %t_pigou(i)=(i-40)/500;
%end
        
for i=1:400
          carbon_price(i)=20*(i-150);
		 %carbon_price(i)=2*(i-150);
		 %t_pigou(i)=(i-40)/500;
end
         
         
for j=1:bs_max
	for inc=1:6
	 if inc==1
	 	set_seed=10;
	 elseif  inc==2
	 	set_seed=3;	
	 elseif  inc==3
	 	set_seed=6;
	 elseif  inc==4
	 	set_seed=4;
	 elseif  inc==5
	 	set_seed=3;
	 elseif  inc==6
	 	set_seed=7;
	 elseif  inc==16
	 	set_seed=1;  
 	 end
 	 
   	%results_file=strcat('\\fsa\faculty\shoude\TempDataCode\EEGap\PolicyAnalysis\Results_PolicyAnalysis\','Optimal_m_adj_t_pigou_heter_em_fac_v4_EEgap_d0.02_44000_inc_',num2str(inc),'_sd',num2str(set_seed),'_bs_',num2str(bs_inc(j,inc)),'_FOX_bs.mat');	    
   	results_file=strcat(path_results{1},'Optimal_m_pos_t_pigou_heter_em_fac_EEgap_d0.02_44000_inc_',num2str(inc),'_sd',num2str(set_seed),'_bs_',num2str(bs_inc(j,inc)),'_FOX_bs.mat');	    
                                         
    opt_inc=load(results_file);
	DSW_heter_inc{j,inc}=opt_inc.Total_SW - cell2mat(opt_inc.Total_base_SW(bs_inc(j,inc)));
    DW_heter_inc{j,inc}=opt_inc.Total_W_avg - cell2mat(opt_inc.Total_base_W_avg(bs_inc(j,inc)));
    Dext_heter_inc{j,inc}=opt_inc.Total_ext_avg - cell2mat(opt_inc.Total_base_ext_avg(bs_inc(j,inc)));
    
	DSW_hom_inc{j,inc}=opt_inc.Total_hom_SW - cell2mat(opt_inc.Total_base_hom_SW(bs_inc(j,inc)));	
    DW_hom_inc{j,inc}=opt_inc.Total_hom_W_avg - cell2mat(opt_inc.Total_base_hom_W_avg(bs_inc(j,inc)));
    Dext_hom_inc{j,inc}=opt_inc.Total_hom_ext_avg - cell2mat(opt_inc.Total_base_hom_ext_avg(bs_inc(j,inc)));
    
    %DSW_heter_inc{j,inc}=opt_inc.Total_SW ;
    %DW_heter_inc{j,inc}=opt_inc.Total_W_avg ;
    %Dext_heter_inc{j,inc}=opt_inc.Total_ext_avg ;
    
	%DSW_hom_inc{j,inc}=opt_inc.Total_hom_SW ;	
    %DW_hom_inc{j,inc}=opt_inc.Total_hom_W_avg ;
    %Dext_hom_inc{j,inc}=opt_inc.Total_hom_ext_avg ;
    
    end
	
% Aggregate across income groups    
   DSW_heter_mat = [DSW_heter_inc{j,1}', DSW_heter_inc{j,2}', DSW_heter_inc{j,3}', ...
              DSW_heter_inc{j,4}', DSW_heter_inc{j,5}', DSW_heter_inc{j,6}']; 	      
   DSW_heter_cty = DSW_heter_mat*share_inc_heter';     
   DSW_heter(j,:) = mean(DSW_heter_cty,2);
          
   DW_heter_mat = [DW_heter_inc{j,1}', DW_heter_inc{j,2}', DW_heter_inc{j,3}', ...
              DW_heter_inc{j,4}', DW_heter_inc{j,5}', DW_heter_inc{j,6}']; 	      
   DW_heter_cty = DW_heter_mat*share_inc_heter';     
   DW_heter(j,:) = mean(DW_heter_cty,2);
   
   Dext_heter_mat = [Dext_heter_inc{j,1}', Dext_heter_inc{j,2}', Dext_heter_inc{j,3}', ...
              Dext_heter_inc{j,4}', Dext_heter_inc{j,5}', Dext_heter_inc{j,6}']; 	      
   Dext_heter_cty = Dext_heter_mat*share_inc_heter';     
   Dext_heter(j,:) = mean(Dext_heter_cty,2);
   
   
   DSW_hom_mat = [DSW_hom_inc{j,1}', DSW_hom_inc{j,2}', DSW_hom_inc{j,3}', ...
              DSW_hom_inc{j,4}', DSW_hom_inc{j,5}', DSW_hom_inc{j,6}']; 	      
   DSW_hom_cty = DSW_hom_mat*share_inc_heter';     
   DSW_hom(j,:) = mean(DSW_hom_cty,2);
          
   DW_hom_mat = [DW_hom_inc{j,1}', DW_hom_inc{j,2}', DW_hom_inc{j,3}', ...
              DW_hom_inc{j,4}', DW_hom_inc{j,5}', DW_hom_inc{j,6}']; 	      
   DW_hom_cty = DW_hom_mat*share_inc_heter';     
   DW_hom(j,:) = mean(DW_hom_cty,2);
   
   Dext_hom_mat = [Dext_hom_inc{j,1}', Dext_hom_inc{j,2}', Dext_hom_inc{j,3}', ...
              Dext_hom_inc{j,4}', Dext_hom_inc{j,5}', Dext_hom_inc{j,6}']; 	      
   Dext_hom_cty = Dext_hom_mat*share_inc_heter';     
   Dext_hom(j,:) = mean(Dext_hom_cty,2);
   
%Identify the max          
   [DSW_heter_max(j), max_id_heter] = max(DSW_heter(j,:));
	carbon_price_star_heter(j) = carbon_price(max_id_heter);       
    DW_heter_max(j) = DW_heter(j,max_id_heter);   
    Dext_heter_max(j) = Dext_heter(j,max_id_heter);       
          
   [DSW_hom_max(j), max_id_hom] = max(DSW_hom(j,:));
	carbon_price_star_hom(j) = carbon_price(max_id_hom);       
    DW_hom_max(j) = DW_hom(j,max_id_hom);   
    Dext_hom_max(j) = Dext_hom(j,max_id_hom);
        
    
end

%Take the mean and s.e. across bootstrap iterations
carbon_price_star_heter_mean = sum(carbon_price_star_heter)/(bs_max);
carbon_price_star_heter_se = std(carbon_price_star_heter,1)/((bs_max-1)^.5);

carbon_price_star_hom_mean = sum(carbon_price_star_hom)/(bs_max);
carbon_price_star_hom_se = std(carbon_price_star_hom,1)/((bs_max-1)^.5);

DSW_heter_mean  = sum(DSW_heter_max)/(bs_max);
DSW_heter_se    = std(DSW_heter_max,1)/((bs_max-1)^.5);

DW_heter_mean  = sum(DW_heter_max)/(bs_max);
DW_heter_se    = std(DW_heter_max,1)/((bs_max-1)^.5);

Dext_heter_mean  = sum(Dext_heter_max)/(bs_max);
Dext_heter_se    = std(Dext_heter_max,1)/((bs_max-1)^.5);

DSW_hom_mean  = sum(DSW_hom_max)/(bs_max);
DSW_hom_se    = std(DSW_hom_max,1)/((bs_max-1)^.5);

DW_hom_mean  = sum(DW_hom_max)/(bs_max);
DW_hom_se    = std(DW_hom_max,1)/((bs_max-1)^.5);

Dext_hom_mean  = sum(Dext_hom_max)/(bs_max);
Dext_hom_se    = std(Dext_hom_max,1)/((bs_max-1)^.5);

Results_table=[carbon_price_star_heter_mean, DSW_heter_mean, DW_heter_mean, Dext_heter_mean;...
			   carbon_price_star_heter_se, DSW_heter_se, DW_heter_se, Dext_heter_se;...
			   carbon_price_star_hom_mean, DSW_hom_mean, DW_hom_mean, Dext_hom_mean;...
			   carbon_price_star_hom_se, DSW_hom_se, DW_hom_se, Dext_hom_se];


end


function combs = nmultichoosek(values, k)
%// Return number of multisubsets or actual multisubsets.
if numel(values)==1 
    n = values;
    combs = nchoosek(n+k-1,k);
else
    n = numel(values);
    combs = bsxfun(@minus, nchoosek(1:n+k-1,k), 0:k-1);
    combs = reshape(values(combs),[],k);
end
end



